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Abstract 

We present a new numerical scheme for one dimensional dynamical 
systems. This is a modification of the discrete gradient method and 
keeps its advantages, including the stability and the conservation of 
the energy integral. However, its accuracy is higher by several orders 
of magnitude. 
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The discrete gradient methods were introduced many years ago in order 
to integrate numerically iV-body systems of classical mechanics with possible 
applications in molecular dynamics and celestial mechanics pQ (see also [2]). 
In this paper we consider the one-dimensional case: 

p = —V'(x) , x — p , (1) 

where V(x) is a potential, and dot and prime denote differentiation with 
respect to t and x, respectively. In this case discrete gradient methods reduce 
to the so called modified midpoint rule: 

Xn+l X n 1 . » 
= g (P»H-1 + Pn) ■ 

(2) 

p n+ l - p n _ V(x n +l) - V(x n ) 
£ Xn+l X n 

One can easily prove that the total energy is preserved by this scheme: 

X -p 2 n + V{x n ) = E = const , (3) 

The modified midpoint rule has been extended, in a natural way, on the 
three-dimensional case and on systems of particles, exactly preserving the 
total energy, the total linear momentum, and the total angular momentum 
of the system pQ. 

More recently discrete gradient methods have been developed in the con- 
text of geometric numerical integration [3], see [HE]. In particular, Quispel 
and his collaborators constructed numerical integrators preserving integrals 
of motion of a given system of ordinary differential equations [HI El E] • 

In general, geometric numerical integrators are very good in preserving 
qualitative features of simulated differential equations but it is not easy to 
enhance their accuracy. Our research is concentrated on improving the effi- 
ciency of the discrete gradient method without loosing its outstanding qual- 
itative advantages. 



In a previous paper we compared several discretizations of the simple 
pendulum equation with a special stress on the long-time behaviour [9] . The 
discrete gradient scheme was among the best ones, especially when large 
energies (rotational motion) and the neighbourhood of the separatrix were 
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concerned. In the paper [H] we proposed a modification of the discrete gra- 
dient scheme ([2]). Assuming the stable equilibrium at x = 0, we replaced e 
by the function 80 = 80(e) 

. 2 u e , s 

8 = — tan — , (4) 

where uq = \/V"(0). The motivation was to preserve (almost exactly) small 
oscillations around x = 0, when the pendulum can be treated as a harmonic 
oscillator. The classical harmonic oscillator admits the so called exact (or 
the best) discretization [TOJ, ITT] . The exact discretization of the harmonic 
oscillator equation was earlier used to derive an orbit-preserving discretiza- 
tion of the classical Kepler problem [12]. The result obtained in the paper 
[5] has been quite satisfying: in the case of small oscillations our method was 
better by 4 orders of magnitude than all other considered schemes (including 
the discrete gradient method). In the case of other initial conditions the 
new method was comparable with the discrete gradient method (by the way: 
both methods are of second order). 

As the main result of this paper we propose another, much more efficient, 
modification of the discrete gradient scheme: 

— 2 = 2 ( Pn+1 + Pn > • 

(5) 

Pn+1 ~ Pn = V(x n+ i) ~ V(x n ) 
8n 2<ri 

where 8 n is a function defined by 

8 n = — tan^, (if V"(x n ) > 0) , 

LO n 2 



J n j 



(if V"(x n ) = 0) , (6) 



8 n = — tanh^, (if V"(x n ) < 0) , 
e denotes the time step, and, finally, 



U n = V\V"(Xn)\ . (7) 
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In fact we have replaced e, appearing in formulas d2J), by a variable 5 n de- 
pending not only on e but also on x n . The crucial point is to choose the best 
form of the function 8 n . In order to obtain ([6]) we have to require that the 
modified scheme (jSJ) is locally exact |13j . It means that the linearization of 
the scheme (|5]) around any fixed x, i.e., 

Cn+l — £n 1 / x 

S " 2 , (8) 

p " + ' " n = -v(x) - i(e„ +1 + f„)v"(x) , 

(where £ n = x n —x) coincides with the exact discretization of the linearization 
of the considered system ([I]) 

! = -"*<*) -V<*){, |=P, 0) 

where x is treated as a constant (equations ([HD, equivalent to the harmonic 
oscillator with a constant force, admit the exact discretization). Thus we 
can express S n in terms of V"(x). Then, identifying x with x n , we arrive at 
formulas (J6]). The detailed derivation of the numerical scheme ([5]) will be 
presented elsewhere [13]. We point out that S n is a function which is almost 
constant (especially for small e), namely: 5 n ~ e, see Fig. [TJ 

The replacement e — > <5 n works very well for the first order system (jSj). 
However, if we try to replace e by 5 n in the second order discrete equation 
for Eq. (46) in the paper [9]), then we get a difference scheme which 

is only a little bit better than the scheme with constant 5 = 5q. 

In practical implementation we use the implicit scheme (jSJ) as the correc- 
tor, while taking as the predictor the explicit scheme 

sm(u n e) 1 — cos(o; n 5) Tr/ . . ... Tjr „, . . 
x n+1 = x n + V > n V\x n ) , (if y"(x„) > 0) , 

x n+1 = x n + ep n - l -e 2 V'{x n ) , (if V"(x n ) = 0) , (10) 

sinh(u; n £) 1 — cosh(a; n £) . ... . 

x n+1 = x n + V\x n ) , (if 7" x B < 0) . 
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In order to obtain (1101) one has to eliminate p n +i from the system §5§ and 
expand the result in the Taylor series with respect to x n+ i — x n , leaving only 
linear terms. 

The numerical scheme fflQj) has the same order (third) as (j3J). However, 
this is of some advantage only for very small e and for very short times (thus 
( TTUj) can serve as a very good predictor). As far as the long-time behaviour 
is concerned the scheme ffTUl) is not good and yields solutions with wrong 
qualitative behaviour. 

The proposed numerical integrator (jSJ) (we will refer to it as locally exact 
discrete gradient method) has important advantages: 

• exact conservation of the energy integral (i.e., eq. (j3J) holds), 

• higher order (third) as compared with the discrete gradient method, 

• high stability and accuracy, 

• very good long-time behaviour of numerical solutions. 

The accuracy of our new numerical scheme was tested on the case of the 
simple pendulum equation (V(x) = — cos a;). We compared it with the stan- 
dard leap-frog scheme, the discrete gradient method and modified discrete 
gradient method (introduced in [5]). All motions of the simple pendulum 
are periodic, so we focus our attention on the relative error of the period 
of considered discretizations. The details of our numerical experiments are 
explained in [9]. Here we just mention that for simplicity we assume xq = 0. 
In this case p$ = IE + 2. 

Both discrete gradient methods studied in [9] are very stable and our 
new method shares this property as well. The accuracy of the new method is 
surprisingly high, especially for small (but not necessarily very small) time- 
steps. As an example we consider the case e = 0.02 (see Fig. For small 
oscillations (e.g., po = 0.02) the accuracy of our new method is greater by 5 
orders of magnitude as compared to the modified discrete gradient method, 
and 9 (nine!) orders of magnitude better than leap-frog or discrete gradient 
method. Actually, the locally exact discrete gradient method is much more 
accurate than any other considered method for any initial conditions, see 
Fig. [2] (note that the scale on the vertical axis is logarithmic). Usually other 
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methods are worse at least by 4 orders of magnitude, except the discrete 
gradient method which, for large p$, is worse "only" by 2 orders of magnitude 
and this difference seems to diminish for larger energies. 

We present also the time-step dependence of the considered numerical 
schemes for three chosen energies (see Fig. [3], Fig. H] and Fig. [51 note that 
these graphs are semilogarithmic). The locally exact discrete gradient scheme 
is almost always the best with very few exceptions. Namely, in the case of 
small oscillations and larger time-steps the modified discrete gradient scheme 
yields similar results, see Fig. [3j Then, for larger amplitudes and larger time- 
steps the discrete gradient scheme is not much worse, see Fig. [5j Finally, the 
leap-frog is comparable with our new method when p is close to a peculiar 
value po ~ 1.21, see Fig. HI This is the somewhat mysterious "resonance 
value" of po for the leap-frog scheme, already reported in [9] . For this initial 
condition the leap-frog method is exceptionally accurate and, for larger time- 
steps, is comparable even with the locally exact discrete gradient scheme. 

The neighbourhood of the separatrix (p « 2) is most difficult to be 
simulated numerically. The discrete gradient method was relatively good at 
this region, see [9]. The locally exact discrete gradient method is excellent 
even in that case, see Fig. [6] and Fig. We are very close to the separatrix 
( |Po — 2 1 = 1CT 10 ), e is large, and nevertheless our new method simulates very 
accurately the motion of the pendulum. Discrete points x n practically lie on 
the continuous curve of the exact solution. The other two discrete gradient 
methods yield good results (at least qualitatively) while the leap-frog scheme 
fails to reproduce even the qualitative behaviour. 

The method presented in this paper can be extended and generalized on 
multidimensional cases and some other numerical integrators (including the 
implicit midpoint rule) [T3]. The main idea of this novel approach consists 
in modifying numerical integrators by replacing e by appropriate functions 
(e.g., 8 n ), depending on e and independent variables, in order to obtain "lo- 
cally exact" integrators (for instance, integrators which are locally equivalent 
to the exact discrete harmonic oscillator). Note that the time-step of the lo- 
cally exact discrete gradient scheme (JSD is equal to e and is assumed to be 
constant. However, there are no obstacles to use the variable time-step and 
it is worthwhile to examinate this possibility in the future. 
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Figure 1: 8 n /e for e = 0.02 (dashed line) and e = 0.1 (solid line). 
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Figure 2: Relative error of the period as a function of po for e = 0.02. White 
triangles: leap-frog, white diamonds: discrete gradient, black diamonds: modified 
discrete gradient (5 = const), black squares: locally exact discrete gradient. 
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Figure 4: Relative error of the period as a function of e for pq = 1.21 ("resonance 
value" for the leap-frog scheme). Symbols: the same as in figure [2j 
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Figure 6: function of n, very near the separatrix (pq = 1.9999999999), for 

e = 0.9. Symbols: the same as in figure [2j The solid line corresponds to the exact 
(continuous) solution. 
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Figure 7: function of n, very near the separatrix (po = 2.0000000001), for 

e = 0.7. Symbols: the same as in figure [2J The solid line corresponds to the exact 
(continuous) solution. 
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